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Abstract 



We present a computer program for the simulation of Mie scattering in case of 
arbitrarily large size parameters. The elements of the scattering matrix, efficiency 
factors as well as the corresponding cross sections, the albedo and the scattering 
asymmetry parameter are calculated. Single particles as well as particle ensembles 
consisting of several components and particle size distributions can be considered. 

O ' PACS: 94.10.Gb Absorption and scattering of radiation 

c/3 . Key words: Mie scattering, scattering matrix, efficiency factors, cross sections, 



radiative transfer - astrophysics, optics, geophysics, biophysics 
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1 Program summary 

Title of program: miex 

Catalogue identifier: 

Program obtainable from: CPC Program Library, Queen's University of Bel- 
fast, N. Ireland 

Computer for which the program is designed and others on which it has been 
tested: 

Computers: Any machine running standard FORTRAN 90; miex has been 
tested on an Intel Celeron processor (Redhat Linux 9.0, Intel Fortran Compiler 
7.1), an Intel XEON processor (SuSE Linux 9.0, Intel Fortran Compiler 8.0), 
and a Sun-Blade-1000 (OS 8.5, Sun Workshop Compiler Fortran 90 2.0). 

Installations: standard 

Operating systems or monitors under which the program has been tested: 
Redhat Linux 9.0, SuSE Linux 9.0, Sun OS 8.5 

Programming language used: Fortran 90 

Memory required to execute with typical data: 1 MByte - several 100 MByte 
(see Sect. A for examples) 

No. of bits in a word: 8 

No. of processors used: 1 

Has the code been vectorized or parallelized ?: NO 

No. of bytes in distributed program, including test data, etc.: 

Distributed format: ASCII 

Keywords : Mie scattering, size parameter, particle ensemble, scattering ma- 
trix, Stokes vector, efficiency factor, cross section, asymmetry parameter, albedo 

Nature of the physical problem: Among a variety of applications, Mie scatter- 
ing is of essential importance for the continuum radiative transfer in cosmic 
dust configurations. In this particular case, Mie theory describes the inter- 
action of electromagnetic radiation with spherical dust grains on the basis 
of their complex refractive index and size parameter. Both, broad grain size 
distributions (radii a: nanometers - millimeters) and a very wide wavelength 



range (A ~ 10~ 10 — 10~ 2 m) of the interacting radiation are considered. Pre- 
vious numerical solutions to the Mie scattering problem are not appropriate 
to consider size parameters x = 2ira/\ > 10 4 — 10 5 . In contrast to this, the 
presented code allows to consider arbitrary size parameters. It will be useful 
not only for applications in astrophysics but also in other fields of science 
(atmospheric and ocean optics, biophysics, etc.) and industry (particle sizing, 
ecology control measurements, etc.). 

Method of solution: Calculations of Mie scattering coefficients and efficiency 
factors as outlined by Voshchinnikov (2004), combined with standard solutions 
of the scattering amplitude functions. Single scattering by particle ensembles 
is calculated by proper averaging of the respective parameters. 

Restrictions on the complexity of the problem: Single Scattering 

Typical running time: Seconds to minutes. 

Unusual features of the program: None 

References: 

Bohren C.F., Huffman D.R., Absorption and scattering of light by small par- 
ticles. John Wiley & Sons, New York (1983) 

Voshchinnikov N.V.: "Optics of Cosmic Dust. I", Astrophysics and Space 
Physics Reviews 12, 1 (2004) 



2 Introduction 



The motivation for the development of the presented Mie scattering routine 
was given by the goal to model the continuum radiative transfer in the circum- 
stellar environment of young stellar objects and X-ray halos where dust scat- 
tering and thermal dust reemission dominate the continuum radiation field. 
In the most approaches, dust grains are assumed to be homogeneous spheres. 
Their optical properties are defined through Mie theory with following input 
parameters: (a) the complex refractive index m = n + ik of a material relative 
to the surrounding medium and (b) the dimensionless size parameter 

x = ^, (1) 



connecting the particle radius a and the wavelength of the incident electro- 
magnetic radiation A. New findings suggest that the radii of the largest grains 
in circumstellar disks around young stellar objects may reach several millime- 
ters to centimeters in radius (see, e.g., Calvet et al. 2002). Since the embedded 
stars emit a considerable amount of energy even in the ultraviolet wavelength 
range (effective temperature of a typical T Tauri star: T c s =4000 K - Gullbring 
et al. 1998), the resulting size parameter x amounts up to 10 5 — 10 6 . In more 
evolved circumstellar disks, so-called "debris disks" , the size parameter x may 
even exceed this value by several orders of magnitude. Furthermore, interstel- 
lar dust grains located near the line of sight toward distant X-ray sources will 
scatter some of its radiation and thus produce a diffuse halo around a point- 
like source. Observations are performed at energies E ~ lOkeV (or more see, 
e.g., Woo et al., 1994) that correspond to wavelengths A ~ 1.24 A. In this 
domain, the size parameter is equal to x = 5067(a/fim)(E/keV) and exceeds 
5 x 10 6 if a > 100 /im. 

For the solution of the described radiative transfer problems, the efficiency 
factors and scattering matrix are required. Under the assumption of spheri- 
cal, homogeneous particles, these quantities are calculated from the Mie series 
solution. The most widely applied numerical Mie code is that from Bohren & 
Huffman (1983). It is based on the improved algorithm of Mie scattering of 
Wiscombe (1980, 1996) and is restricted to size parameters < 2 x 10 4 . There- 
fore, this numerical solution is not appropriate for the applications described 
above (see also Shah 1977, 1992 for Mie calculations in the size parameter 
regime x < 10 5 ). 

The computer program presented in this article is based on the calculation 
of the Mie coefficients and efficiency factors for single particles as described 
by Voshchinnikov (2004). The calculations of the angular functions are based 
on the standard approach (see Wiscombe 1980, Bohren & Huffman 1983). We 



successfully tested the direct implementation of Legendre functions (Zhang & 
Jin 1996) for the determination of the scattering amplitude functions as well, 
but only the fastest, i.e., the first approach is realized in the published version 
of the code. To account for particle size distributions and/or ensembles of 
different components, a proper averaging of the single parameters (efficiency 
factors etc.) is implemented as well. 



3 Long Write-up 



3.1 Mie scattering 



In the following, a brief overview of the scattering process and the different 
quantities which can be derived with miex is given. For a more detailed de- 
scription of the Mie scattering solution we refer to the books of van de Hulst 
(1957) and Bohren & Huffman (1983). 

Let the radiation be described by the four-component Stokes vector I = 
(J, Q, U, V) T . The scattering of radiation by a spherical particle is described 
by the scattering (Miiller) matrix F of special type 
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where Iq (ii) is the Stokes vector before (after) the scattering and 
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Here, is the angle between the direction of the incident and the scattered 
radiation (scattering angle). The elements of the scattering matrix i 7 ^ can 
be derived from the complex amplitude functions (the asterisk denotes the 
complex conjugation) Si, S 2 : 
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F 12 (e) = ±[\s 2 (e)\ 2 -\Sim 



F n (e) = -[\S 2 (Q)\ Z + \Si(&)\ 

^33(6) = \ [S* 2 (e)S 1 (Q) + S 2 (Q)Sl(Q)} 



(4) 



^34(6) = - [Si(e)s* 2 (e) - s 2 (e)sj(e)] 



The amplitude functions can be calculated as follows: 



°° 277 + 1 

f^i n(n + 1) 



(5) 



and 
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The complex scattering (Mie) coefficients a n and b n depend on the size pa- 
rameter x and the refractive index m of the material. The angular functions 
TT n and T n depend on cos only and can be found from recurrence relations 
(Wiscombe 1980; n > 2) 
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and 



r n (0) = ncos07r„(0) - (n+ l)7T„_i(0). 



The initial values are 



7T O (0) = O, 7ri(6) = l. 



(9) 



The coefficients 
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are usually transformed to a form convenient for calculations. Following Deir- 
mendjian (1969) and Loskutov (1971), we replace the Riccati-Bessel functions 
ip n (x) and ( n (x) and their first derivatives ip' n (x) and C n (x) by spherical Bessel 
functions of the first and second kind J n +i/2(%) and Y n+ \/2{x). With the aid of 
recursive relations for these functions (see, e.g., Abramowitz and Stegun 1964) 
we obtain: 
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The logarithmic derivative to the Bessel functions of a complex argument ^4„ 
is calculated via backward recursion using the relation 
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The starting number (variable "num" in the program) for x < 50000 (A n = 0) 
is chosen according to the recommendation of Loskutov (1971). It is smaller 
for large arguments than that given by Wiscombe (1980). 

According to Voshchinnikov (2004), the calculation of the Bessel functions of 
a real argument is based on the upward recursion for the functions of the 
second kind Y n+ i/ 2 {x), which is known to be stable. This is given by the 
relation (Abramowitz and Stegun 1964; n > 0): 
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with the initial values 
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The Bessel functions of the first kind J n +i/2( x ) (n > 1) are determined by the 
relation 

r 2 1 

J n +i/2(x) = VY n+l/2 (x)J n _ 1/2 (x) Y n _ 1/2 (x), (18) 
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where 
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We want to remark that equation (18) gives wrong results if n ^> x, but this 
is not the case for standard Mie calculations. In order to avoid an overflow in 
the forward recursion, a normalization is used. The description of the general 
method of calculations of Bessel functions can be found in Loskutov (1971). 

The extinction, scattering, backscattering and radiation pressure efficiency 
factors (Qextj Qsca,, Qbk and Q w ) are given by the relations 
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Following quantities can be derived from the above efficiency factors: 



Absorption efficiency : Q abs = Q cxt - Qsca, 

Albedo : A = Q sca /Qcxt, 

Asymmetry parameter : g = (Q cxt - Q pr )/Q sca . 



(24) 



The asymmetry parameter g (or "mean cosine" (cos©)) describes the distri- 
bution of the scattered radiation in the forward / backward direction. It is 
defined as 
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and can be also found from the integration of the element Fn of the scattering 
matrix. 



The corresponding cross sections can be derived from the relation 

C = GQ, (26) 

where G = ira 2 is the geometrical cross section of the particle. 

For the description of scattering by an ensemble of particles consisting of 
several species and particles sizes, weighted mean values of the different quan- 
tities described above can be defined. They are obtained as the sums of the 
corresponding quantity averaged over the corresponding size distribution. As- 
suming "J" different species in the mixture having the fractional abundances 
fj and a particle number density with some size distribution n(a), we can 
formulate the following normalization condition: 



Y,fj I n(a)da = l, ££ = 1, (27) 



where a m i n and a max are the minimum and maximum particle radius. The 
Stokes parameters as well as all cross sections are additive. Therefore, the 
ensemble averaged values can be derived from their weighted contributions 
(see, e.g., Martin 1978, Sole 1980): 
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For the albedo A and the asymmetry parameter g the following expressions 
must be used: 
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3.2 Code description 



The presented Mie scattering code is written in Fortran 90/95. The main pro- 
gram regulates the input / output of data, calls the routine for the calculation 
of Mie scattering characteristics for a single size parameter and particle com- 
position (embedded in the module mie-routines) , and performs the averaging 
in case of a polydisperse ensemble. 

The calculations are performed with double precision accuracy using the in- 
ternally defined data type r2 

integer, parameter, public :: r2 = selected_real_kind(9) (32) 

in the module datatype. Changing the definition of this data type allows one 
to adapt the variables to any required number of significant digits. 

The routines shexqnn2 and aa2 are adaptations of the Mie scattering code 
published by Voshchinnikov (2004). They were extended in order to calculate 
the amplitude functions (see Eqs. (5), (6)) and to satisfy Fortran 90/95 stan- 
dards. The numerical realization for calculation of the scattering amplitude 
functions follows the standard approach (Eqs. (7)-(9)). We also successfully 
tested the implementation of a direct calculation of the Legendre functions, 
based on the routine mlpmn.f or provided by Zhang & Jin (1996), but decided 
not to include it in the presented version of the code because it resulted in 
longer runtimes. The runtime of the code is inversely proportional to the step 
of the scattering angles for which the amplitude functions and, based on 
them, the scattering matrix elements are calculated. The angular step size AG 
is an input parameter. 

The calculation of Mie series for a single particle radius for a given material is 
finished as soon as the relative contribution of the current term to the extinc- 
tion efficiency becomes smaller than 10~ 15 . If the default maximum number 
of terms (2 x 10 7 ) is too small to achieve this accuracy, it may be increased in 
subroutine shexqnn2: 

! Maximum number of terms to be considered 
nterms = 20 000 000 

In order to derive the optical properties of the particle ensemble weighted over 
a size distribution, an arbitrary number of particle radii to be considered can 
be defined. The considered radii are equidistantly distributed on a logarithmic 
scale within the radius interval [o m i n , a max ]. In the published version, the size 
distribution follows a power-law 

n(a) oc a q , (33) 
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where n(a) is the relative number of particles with the radius a and q is a con- 
stant, usually negative quantity for a given component of the mixture. This 
size distribution was introduced by Mathis et al. (1977) for silicate-graphite 
mixture with q = —3.5 in the process of the interpretation of the interstellar 
extinction curve and has been frequently used in different astrophysical appli- 
cations. The code may be easily adapted to other particle size distributions by 
modification of the following program line (representing the Eq. (27) combined 
with Eq. (33)) 

weight = abun(jcomp) * rad * *q * delrad, (34) 



where weight is the weight for the component ^jcomp with the radius rad 
and relative abundance abun. The quantity delrad is the radius step width 
at the current particle radius rad during the numerical integration stated in 
Eq. (27). 

An arbitrary number of chemically different components in the ensemble can 
be considered. Beside the optical properties, each component is characterized 
by its relative abundance in the particle ensemble. The complex refractive 
index as a function of wavelength has to be provided in tabular form for each 
component in a separate file in the directory ./ri-data/ (see Tabl. 1): 

Column 1 : Wavelength A [//m], 
Column 2 : Re{m(A)}, 
Column 3 : Im{m(A)} 

The results are stored in the directory ./results/ (see Tabl. 1, 2). 

The order of the scattering matrix elements in the respective files (see Table 2) 
is defined by the following algorithm: 

do for all wavelengths A [fJ.ni] , in increasing order 

write A 

do for all angles 0; in increasing order 
write [°], F lk (Q,\) 

end do 
end do 
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Table 1 

Directory Structure. 



Directory 



Contents 



•/ 
./ri-data 

./results 



Executable (miex), Source code, Makefile 
Tables with optical data 
Output files (see Tab. 2) 



Table 2 

Output files (project is the project name) 



Filename 



Contents 



project 

project . qext 
project, cext 
project, qsca 
project, csca 
project . qbk 
project . cbk 
project, qabs 
project, cabs 
project . alb 
project, g 
project . qpr 
project . f 1 1 
project, f 12 
project. f 33 
project, f 34 



Main file containing all results 
(see Appendix A for an example) 



A [/mi] 


Vcxt 


A [/.<m] 


Cext [m 2 ] 


A [/.mi] 


Vsca 


A [/mi] 


Csca [m 2 ] 


A [/mi] 


Qbk 


A [/mi] 


C bk [m 2 ] 


A [/mi] 


Vabs 


A [/mi] 


Cabs [m 2 ] 


A [/mi] 


albedo 


A [/mi] 


asymmetry parameter g 


A [/./m] 


Vpr 


©r],Fn(A[/mi],©[ ]) 


©[°],F 12 (A[/mi],©[°]) 


e[°],F 3 3(A[/mi],0[ o ]) 


©[ ],F 3 4(A[/mi],©[ ]) 
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A Test Run 



In the following, miex is used to calculate scattering matrix elements, effi- 
ciency factors, cross sections, the albedo and the asymmetry parameter for 
a particle ensemble typical for astrophysical applications (see Table A.l and 

A.2): 

(1) Three components: 

(a) Astronomical silicate (relative abundance: 62.5%; 
input file: silicate) 

(b) Graphite; refractive index for the electric field vector parallel to the 
crystallographic c axis (relative abundance: 25%; input file: 
grap-per) 

(c) Graphite; refractive index for the electric field vector perpendicular 
to the crystallographic c axis (relative abundance: 12.5%, input file: 
grap-par) 

(2) Grain size distribution: 

• minimum/maximum grain size: a m i n =0.005//m, a max =100/im 

• n(a) oc a -3 ' 5 

• number of considered, logarithmically equidistantly distributed radii: 
100 

(3) Calculation of the scattering matrix elements for the angles 0°, 2°, 4°, 
. . . , 180° 

Input files (with the dust properties) must be located in the directory " . /input/" 
while the results ared stored in the directory "./results/". 

In a second example, documented in Table A. 3 and A. 4, a dust grain ensemble 
with the chemical components as before but with a single grain size of 10 cm 
is considered. 

The calculations are carried out for 100 wavelengths distributed equidistantly 
on a logarithmic scale within the interval [0.05/im, 2000/im]. The optical con- 
stants have been interpolated at the particular wavelengths using the data 
from Weingartner & Draine (2001) and Laor & Draine (1993). 

We want to remark, that a large variety of refractive indexes can be found in 
the Jena - Petersburg Database of Optical Constants (Henning et al. 1999, see 
also Voshchinnikov 2004). The provided dust parameter files are presented in 
the format required by miex. 
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Table A.l 

Input dialog for "examplel". The runtime of the code on an Intel XEON CPU 
3.0 GHz, using the Intel Fortran Compiler 8.0, amounts to 19 seconds (required 
memory: 1.5 MByte). 

Real refractive index of the surrounding medium : 1.0 

Number of wavelengths : 100 

Number of components : 3 

Name of the dust data files (lambda/n/k data) 

001. component : silicate 

002. component : grap-per 

003. component : grap-par 

Relative abundances of the different components [%] 

001. component : 62.5 

002. component : 25 

003. component : 12.5 
-1- Single grain size 

-2- Grain size distribution : 2 

Minimum grain size [micron] : 0.005 

Maximum grain size [micron] : 100.0 

Size distribution exponent : -3.5 

Number of size bins : 100 

Calculate scattering matrix elements (0=n/l=y) : 1 
Number of scattering angles in the interval 
[0°,180°]; odd number! 

[example: '181' -> step width = 1°] : 91 

Project name (8 characters) : examplel 

Save results in separate files (0=n/l=y) : 1 
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Table A.2 

Program output. File: ' ' . /results/examplel' ' 



# *** PROJECT PARAMETERS *** 

# Number of wavelengths : 100 



# *** RESULTS *** 

# 1. Wavelength [micron], Extinction efficiency factor / cross section [m**2] 
5.000000000000000E-002 1.32305406861068 4.681097343048097E-016 

# 2. Wavelength [micron], Scattering efficiency factor / cross section [m**2] 
5.000000000000000E-002 0.587523440006963 2.078716568917760E-016 

# 3. Wavelength [micron], Backscattering efficiency factor / cross section [m**2] 
5.000000000000000E-002 6.859423237184992E-002 2.426932402252230E-017 

# 4. Wavelength [micron], Absorption efficiency factor / cross section [m**2] 
5.000000000000000E-002 0.735530628603716 2.602380774130337E-016 

# 5. Wavelength [micron], Albedo 
5.000000000000000E-002 0.444066084633951 

# 6. Wavelength [micron], Scattering asymmetry factor g 
5.000000000000000E-002 0.800843488171790 

# 7. Radiation pressure efficiency factor Qpr 
5.000000000000000E-002 0.852539747532814 

# 8. Wavelength [micron], theta [degree], F11-F12-F33-F34 
5.000000000000000E-002 0.000000000000000E+000 105571.611539073 
5.000000000000000E-002 0.000000000000000E+000 -5.641242912656195E-012 
5.000000000000000E-002 0.000000000000000E+000 105571.611539073 
5.000000000000000E-002 0.000000000000000E+000 7.328649202947496E-011 



Table A.3 

Input dialog for "example2". The runtime of the code on an Intel XEON CPU 
3.0 GHz, using the Intel Fortran Compiler 8.0, amounts to 129 seconds (required 
memory: ~240 MByte). 

Real refractive index of the surrounding medium : 1.0 

Number of wavelengths : 100 

Number of components : 3 

Name of the dust data files (lambda/n/k data) 

001. component : silicate 

002. component : grap-per 

003. component : grap-par 

Relative abundances of the different components [%] 

001. component : 62.5 

002. component : 25 

003. component : 12.5 
-1- Single grain size 

-2- Grain size distribution : 1 

Grain size [micron] : 100000 

Calculate scattering matrix elements (0=n/l=y) : 

Project name (8 characters) : example2 

Save results in separate files (0=n/l=y) : 1 
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Table A.4 

Program output. File: ' ' . /results/example2' ' 



# *** PROJECT PARAMETERS *** 
■$ Number of wavelengths : 100 



# *** RESULTS *** 

# 1. Wavelength [micron], Extinction efficiency factor / cross section [m**2] 
5.000000000000000E-002 1.34425083271595 4.223088540642379E-002 

# 2. Wavelength [micron], Scattering efficiency factor / cross section [m**2] 
5.000000000000000E-002 0.731763632480501 2.298903251964923E-002 

# 3. Wavelength [micron], Backscattering efficiency factor / cross section [m**2] 
5.000000000000000E-002 0.174029454773743 5.467296566254283E-003 

# 4. Wavelength [micron], Absorption efficiency factor / cross section [m**2] 
5.000000000000000E-002 0.612487200235445 1.924185288677456E-002 

# 5. Wavelength [micron], Albedo 
5.000000000000000E-002 0.544365392731083 

# 6. Wavelength [micron], Scattering asymmetry factor g 
5.000000000000000E-002 0.904868766536477 

^t 7. Radiation pressure efficiency factor Qpr 
5.000000000000000E-002 0.682100777197063 
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